
// C++ ncludes
#include <string>
#include <fstream>
#include <map>
#include <set>

#include "petsc.h"

// Local includes
#include "perf_log.h"
#include "equation_systems.h"
#include "mesh_base.h"
#include "elem.h"

#include "parallel.h"

#include "moab/Core.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/Range.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/CN.hpp"

#include "iMesh.h"
#include "MBiMesh.hpp"

#include "projection_errors.h"
#include "stacktrace.h"
#include "moab_libmesh_wrapper.h"

#define PRINT_LAST_ERROR \
    if (MB_SUCCESS != result) {\
      std::string tmp_str;\
      std::cout << "Failure; message:" << std::endl;\
      mdbImpl->get_last_error(tmp_str);\
      std::cout << tmp_str << std::endl;\
      print_stacktrace(); \
      CHKERRABORT(PETSC_COMM_WORLD, result); }


MOAB_libMesh_Wrapper::MOAB_libMesh_Wrapper (Core* impl, unsigned int num_systems )
        : mdbImpl(impl), n_systems(num_systems)
{
  root_sets.assign(num_systems, 0) ;
  partition_handles.assign(num_systems, 0) ;
  parallel_communicators.resize(num_systems, 0);

  // read it in via the iMesh object
  imesh = reinterpret_cast<iMesh_Instance>( new MBiMesh(mdbImpl) );

  // create the pcomm tag
  Tag pcomm_tag = moab::ParallelComm::pcomm_tag(mdbImpl, true) ;

  int result;
  for (unsigned int isystem=0; isystem < num_systems; ++isystem)
  {
    result = mdbImpl->create_meshset( MESHSET_SET, root_sets[isystem] ); PRINT_LAST_ERROR;
    result = mdbImpl->create_meshset( MESHSET_SET, partition_handles[isystem] ); PRINT_LAST_ERROR;
    parallel_communicators[isystem] = moab::ParallelComm::get_pcomm( mdbImpl, partition_handles[isystem], &PETSC_COMM_WORLD );
  }

  if (!pcomm_tag)
  {
    PetscPrintf(PETSC_COMM_WORLD, "\nParallel Comm Tag doesn't exist. Bug possibly.\n");
    std::cin.get();
  }
//  parallel_communicators[isystem] = moab::ParallelComm::get_pcomm(mdbImpl, partn_handle, &PETSC_COMM_WORLD);
//  parallel_communicators[isystem] = new moab::ParallelComm(mdbImpl, PETSC_COMM_WORLD);
}


MOAB_libMesh_Wrapper::~MOAB_libMesh_Wrapper ()
{
  for (unsigned int i=0; i < n_systems; ++i)
  {
    delete parallel_communicators[i];  parallel_communicators[i] = NULL;
  }
  parallel_communicators.clear();

  mdbImpl->delete_entities(&partition_handles[0], partition_handles.size());
  mdbImpl->delete_entities(&root_sets[0], root_sets.size());

  partition_handles.clear();
  root_sets.clear();
  n_systems = 0;
  mdbImpl = NULL;
}


void MOAB_libMesh_Wrapper::cell_connectivity (const Elem* elem, std::vector<unsigned int>& cell_connectivity)
{
  //  std::vector<unsigned int> ordered_connectivity;
    //elem->connectivity(0,VTK,ordered_connectivity);
    const unsigned int n_nodes = elem->n_nodes();
    cell_connectivity.assign(n_nodes, 0);
    if (n_nodes == 8)   // Linear HEX
    {
        for (unsigned int i=0; i < 8; ++i)
          cell_connectivity[i] = elem->node(i);
    }
    else if (n_nodes == 27)   // Quadratic HEX
    {
        const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
                                        15,16,19,17,18,20,21,24,22,23,25,26};
        const unsigned int moab_nodes[] = {0,1,2,3,4,5,6,7, // corners
                                           8,11,12,9,13,10,14,15,16,19,17,18, // edges
                                           24, 20, 23, 21, 22, 25,                         // faces
                                           26 };

        for (unsigned int i=0; i < 27; ++i)
          cell_connectivity[moab_nodes[i]] = elem->node(nodes[i]);
     }
     else
     {
        PetscPrintf(PETSC_COMM_WORLD, "\nUnknown element type. Please define in cell_connectivity.\n");
		print_stacktrace();
		CHKERRABORT(PETSC_COMM_WORLD, 1);
    }
}


EntityType MOAB_libMesh_Wrapper::cell_type(const Elem* elem)
{
  EntityType celltype = MBMAXTYPE ; // initialize to something to avoid compiler warning

  switch(elem->type())
  {
  case EDGE2:           //VTK_LINE
  case EDGE3:           //VTK_QUADRATIC_EDGE
  case EDGE4:           //VTK_HIGHER_ORDER_EDGE
    celltype = MBEDGE;
    break;// 1
  case TRI3:            //VTK_TRIANGLE
  case TRI6:            //VTK_QUADRATIC_TRIANGLE
    celltype = MBTRI;
    break;// 4
  case QUAD4:          //VTK_QUAD
  case QUAD8:          //VTK_QUADRATIC_QUAD
  case QUAD9:          //VTK_BIQUADRATIC_QUAD
    celltype = MBQUAD;
    break;// 6
  case TET4:           //VTK_TETRA
  case TET10:          //VTK_QUADRATIC_TETRA
    celltype = MBTET;
    break;// 9
  case HEX8:            //VTK_HEXAHEDRON
  case HEX20:
  case HEX27:           //VTK_QUADRATIC_HEXAHEDRON
    celltype = MBHEX;
    break;// 11
  case PRISM6:           //VTK_WEDGE
  case PRISM15:
  case PRISM18:          //VTK_HIGHER_ORDER_WEDGE
    celltype = MBPRISM;
    break;// 13
  case PYRAMID5:
    celltype = MBPYRAMID;   //VTK_PYRAMID
    break;// 14
  case INFEDGE2:
  case INFQUAD4:
  case INFQUAD6:
  case INFHEX8:
  case INFHEX16:
  case INFHEX18:
  case INFPRISM6:
  case INFPRISM12:
  case NODEELEM:
  case INVALID_ELEM:
  default:
        PetscPrintf(PETSC_COMM_WORLD, "\nUnknown element type '%D'. Please define in cell_type.\n", elem->type() );
		print_stacktrace();
		CHKERRABORT(PETSC_COMM_WORLD, 1);
  }

  return celltype;
}

EntityHandle MOAB_libMesh_Wrapper::make_set( unsigned int options,
                         Range& entities,
                         bool reverse,
                         const int id )
{
  ErrorCode result;
  EntityHandle handle;
  result = mdbImpl->create_meshset( options, handle );  PRINT_LAST_ERROR;

  if (reverse)
  {
    moab::Range::iterator rit;
    EntityHandle eh;
    for (Range::reverse_iterator rit = entities.rbegin(); rit != entities.rend(); ++rit) {
    //for (rit = entities.end()-1; rit != entities.begin(); rit--) {
      eh = *rit;
      result = mdbImpl->add_entities( handle, &eh, 1 );  PRINT_LAST_ERROR;
    }
    eh = *rit;
    result = mdbImpl->add_entities( handle, &eh, 1 );  PRINT_LAST_ERROR;
  }
  else
  {
    result = mdbImpl->add_entities( handle, entities ) ;  PRINT_LAST_ERROR;
  }

  if (id)
  {
    Tag id_tag;
    result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag,
                                  MB_TAG_CREAT|MB_TAG_DENSE ); PRINT_LAST_ERROR;

    result = mdbImpl->tag_set_data( id_tag, &handle, 1, &id ) ; PRINT_LAST_ERROR;
  }

  return handle;
}

// Create an entity set containing the passed entities.
// options      - set options to pass to create_meshset
// entities     - array of entities to put in set
// num_entities - length of 'entities'
// reverse      - if true, add entities in reverse order
// id           - value for global id on set
EntityHandle MOAB_libMesh_Wrapper::make_set( unsigned int options,
                         EntityHandle* entities,
                         size_t num_entities,
                         bool reverse,
                         const int id )
{
  ErrorCode result;
  EntityHandle handle;
  result = mdbImpl->create_meshset( options, handle );  PRINT_LAST_ERROR;

  if (num_entities > 0)
  {
    if (reverse)
    {
      for (int i = (int)num_entities - 1; i >= 0; --i)
      {
        result = mdbImpl->add_entities( handle, entities + i, 1 );  PRINT_LAST_ERROR;
      }
    }
    else
    {
      result = mdbImpl->add_entities( handle, entities, num_entities ) ;  PRINT_LAST_ERROR;
    }
  }

  if (id)
  {
    Tag id_tag;
    result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, id_tag,
                                  MB_TAG_CREAT|MB_TAG_DENSE ); PRINT_LAST_ERROR;

    result = mdbImpl->tag_set_data( id_tag, &handle, 1, &id ) ; PRINT_LAST_ERROR;
  }

  return handle;
}

ErrorCode MOAB_libMesh_Wrapper::write_equation_systems(const std::string& fname, const EquationSystems& physics, unsigned int isystem,
                            bool build_solutions, bool build_adjacencies, bool verbose)
{
	ErrorCode result = MB_SUCCESS;
	int nprocs, rank;
  std::string writeOpts = "" ;

  assert(isystem < n_systems);

  MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  const MeshBase& mesh = (MeshBase&)physics.get_mesh();
  unsigned int dim = mesh.spatial_dimension();
  MeshBase::const_element_iterator element_it = mesh.active_elements_begin() ;
  const MeshBase::const_element_iterator element_end = mesh.active_elements_end();
  MeshBase::const_node_iterator node_it = mesh.nodes_begin() ;
  const MeshBase::const_node_iterator node_end = mesh.nodes_end();

  const unsigned int n_nodes = mesh.n_nodes();
  const unsigned int n_elems = mesh.n_elem();
  unsigned int index = 0;

  std::vector<EntityHandle> global_nodes, global_elems;

  //parallel_communicators[isystem] = moab::ParallelComm::get_pcomm(mdbImpl, partition_handles[isystem], &PETSC_COMM_WORLD);
  //parallel_communicators[isystem] = new moab::ParallelComm(mdbImpl, PETSC_COMM_WORLD);

  //result = mdbImpl->create_meshset( MESHSET_SET, root_sets[isystem] ); PRINT_LAST_ERROR;
  //root_sets[isystem] = make_set( MESHSET_SET, 0, 0, true,  EMPTY_SET_ID*(isystem+1)  );

  if (!rank)
  {
    // Initialize the moab object
    result = mdbImpl->set_dimension(dim);

    global_nodes.assign(n_nodes, 0);
    std::vector<std::vector<EntityHandle> > local_nodes(nprocs);

    // write mesh nodes onto MOAB instance
    double vcoords[3];
    for(node_it = mesh.nodes_begin(); node_it!=node_end; ++node_it)
    {
      const Node& nd = *(*node_it);

      vcoords[0] = nd(0); vcoords[1] = nd(1); vcoords[2] = nd(2);
      result = mdbImpl->create_vertex(vcoords, global_nodes[nd.id()]); PRINT_LAST_ERROR;

      // store the vertex corresponding to a certain partition
      local_nodes[nd.processor_id()].push_back(global_nodes[nd.id()]);
    }
    result = mdbImpl->add_entities ( root_sets[isystem], &global_nodes[0], global_nodes.size() ); PRINT_LAST_ERROR;

    std::vector<unsigned int> cell_connectivity;
    std::vector<EntityHandle> moab_connectivity;

    std::vector<std::vector<EntityHandle> > local_elems(nprocs);
    global_elems.assign(n_elems, 0);

    for(element_it = mesh.elements_begin(); element_it!=element_end; ++element_it)
    {
      const Elem* elem = *element_it;

      this->cell_connectivity(elem, cell_connectivity);
      moab_connectivity.resize(cell_connectivity.size());

      //PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tConnectivity Element %D", rank, index) ;
      for (unsigned int i=0; i < cell_connectivity.size(); ++i) {
        moab_connectivity[i] = global_nodes[cell_connectivity[i]];

      // libMesh connectivity rearrange: HEX27: 0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,15,16,19,17,18,20,21,24,22,23,25,26
      // MOAB connectivity rearrange
/*
              const elementDefinition& eletype = eletypes_imp[type];
              nnodes = eletype.nnodes;

      int *reorder = GmshUtil::gmshElemTypes[curr_elem_type][4];

      // Copy passed element connectivity into entity sequence data.
      if (type.node_order)
      {
        for (unsigned long i = 0; i < num_elem; ++i)
          for (int j = 0; j < node_per_elem; ++j)
            conn_array[i*node_per_elem+type.node_order[j]] = connectivity[i*node_per_elem+j];
      }
      else
      {
        memcpy( conn_array, &connectivity[0], connectivity.size() * sizeof(EntityHandle) );
      }
*/

        //PetscPrintf(PETSC_COMM_SELF, "\t%D %D", cell_connectivity[i], moab_connectivity[i]) ;
      }

//      EntityHandle elem_handle ;
      result = mdbImpl->create_element (cell_type(elem), &moab_connectivity[0], moab_connectivity.size(), global_elems[elem->id()]);PRINT_LAST_ERROR;

      local_elems[elem->processor_id()].push_back(global_elems[elem->id()]);
    }
    result = mdbImpl->add_entities ( root_sets[isystem], &global_elems[0], global_elems.size() ); PRINT_LAST_ERROR;

    Range dim3, dim2, adj;
    if (build_adjacencies)
    {
      PetscPrintf(PETSC_COMM_SELF, "[%D]\tGetting adjacencies.\n", rank);
      mdbImpl->get_entities_by_dimension( root_sets[isystem], 3, dim3, true );

      // generate all lower dimensional adjacencies
      mdbImpl->get_entities_by_dimension( root_sets[isystem], 2, dim2, true );
      mdbImpl->get_adjacencies( dim3, 1, true, adj, Interface::UNION );
      mdbImpl->get_adjacencies( dim2, 1, true, adj, Interface::UNION );
      mdbImpl->get_adjacencies( dim3, 2, true, adj, Interface::UNION );

      mdbImpl->add_entities( root_sets[isystem], adj );
    }

    // create some sets
    EntityHandle vertex_set, region_set, face_set;
    //vertex_set = make_set( MESHSET_ORDERED|MESHSET_TRACK_OWNER, &nodes_vtk_id[node_offset], n_local_nodes, true,  VERTEX_SET_ID  );
  	vertex_set = make_set( MESHSET_ORDERED, &global_nodes[0], global_nodes.size(), false,  MOAB_libMesh_Wrapper::VERTEX_SET_ID*(isystem+1)  );
    if (build_adjacencies)
      face_set   = make_set( MESHSET_SET,     dim2, false, MOAB_libMesh_Wrapper::FACE_SET_ID*(isystem+1)    );
  	region_set = make_set( MESHSET_SET,     &global_elems[0], global_elems.size(), false, MOAB_libMesh_Wrapper::REGION_SET_ID*(isystem+1)  );

    Tag partn_tag = parallel_communicators[isystem]->partition_tag	();

    // next get the global ID tag handle
    Tag dimId = mdbImpl->geom_dimension_tag();

    // next get the global ID tag handle
    Tag globalId = mdbImpl->globalId_tag();
  	result = mdbImpl->tag_set_data(globalId, &global_nodes[0], global_nodes.size(), &global_nodes[0]);PRINT_LAST_ERROR;
    result = mdbImpl->tag_set_data(globalId, &global_elems[0], global_elems.size(), &global_elems[0]);PRINT_LAST_ERROR;

    // set the partition to the global vertex/region set => 1 proc file
    //result = mdbImpl->tag_set_data(partn_tag, &vertex_set, 1, &rank);PRINT_LAST_ERROR;
    //result = mdbImpl->tag_set_data(partn_tag, &region_set, 1, &rank);PRINT_LAST_ERROR;

    std::vector<EntityHandle> vsets(nprocs);
    std::vector<EntityHandle> esets(nprocs);
    for (int irank=0; irank < nprocs; irank++)
    {
      vsets[irank] = make_set( MESHSET_ORDERED, &local_nodes[irank][0], local_nodes[irank].size(), false/*, VERTEX_SET_ID*/  );
      esets[irank] = make_set( MESHSET_ORDERED, &local_elems[irank][0], local_elems[irank].size(), false/*, REGION_SET_ID*/ );

      int local_set_id = MOAB_libMesh_Wrapper::SET_SET_ID*(irank+1)*(isystem+1);
      EntityHandle sets[] = {vsets[irank], esets[irank]};
      make_set( MESHSET_ORDERED, sets,      2, false, local_set_id     );

      //const unsigned empty_flags = MESHSET_ORDERED|MESHSET_TRACK_OWNER;
      //empty_set  = make_set( empty_flags,     0,         0, true,  EMPTY_SET_ID   );
      //EntityHandle sets[] = {face_set, vertex_set, region_set, empty_set};
      //EntityHandle sets[] = { vertex_set, region_set, empty_set};
      //make_set( MESHSET_ORDERED, sets,      3, false, SET_SET_ID     );

      result = mdbImpl->tag_set_data(partn_tag, sets, 2, &irank);PRINT_LAST_ERROR;
      result = mdbImpl->add_entities( root_sets[isystem], sets, 2 ) ;PRINT_LAST_ERROR;

      //result = mdbImpl->tag_set_data(partn_tag, &vsets[irank], 1, &irank);PRINT_LAST_ERROR;
      //result = mdbImpl->tag_set_data(partn_tag, &esets[irank], 1, &irank);PRINT_LAST_ERROR;

      //result = mdbImpl->add_entities( root_sets[isystem], &vsets[irank], 1 ) ;PRINT_LAST_ERROR;
      //result = mdbImpl->add_entities( root_sets[isystem], &esets[irank], 1 ) ;PRINT_LAST_ERROR;
    }
    result = parallel_communicators[isystem]->set_partitioning(root_sets[isystem]);

    int vdim = 0, fdim=dim-1;
    result = mdbImpl->tag_set_data(dimId, &region_set, 1, &dim);PRINT_LAST_ERROR;
    result = mdbImpl->tag_set_data(dimId, &vertex_set, 1, &vdim);PRINT_LAST_ERROR;
  	if (build_adjacencies) {
  	  result = mdbImpl->add_entities( root_sets[isystem], &face_set, 1 ) ;PRINT_LAST_ERROR;
  	  result = mdbImpl->tag_set_data(dimId, &face_set, 1, &fdim);PRINT_LAST_ERROR;
  	}

    if (verbose) {
      moab::Range ents;
      // iterate over dimensions
      for (unsigned int d = 0; d <= dim; d++) {
        ents.clear();
        result = mdbImpl->get_entities_by_dimension((EntityHandle)root_sets[isystem], d, ents);PRINT_LAST_ERROR;
        PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tFound %d %d-dimensional entities\n", rank, ents.size(), d);
        for (moab::Range::iterator it = ents.begin(); it != ents.end(); it++) {
          moab::EntityHandle ent = *it;
          PetscPrintf(PETSC_COMM_SELF, "[%D]\t\tFound d=%d entity %d.\n", rank, d, mdbImpl->id_from_handle(ent));
        }
      }
    }

    result = mdbImpl->add_entities( root_sets[isystem], &region_set, 1 ) ;PRINT_LAST_ERROR;
  	result = mdbImpl->add_entities( root_sets[isystem], &vertex_set, 1 ) ;PRINT_LAST_ERROR;

    // finally add the file set to the mesh's root set
  	//result = mdbImpl->add_entities( 0, &root_sets[isystem], 1 ) ;PRINT_LAST_ERROR;
  	//result = mdbImpl->add_parent_meshset( root_sets[isystem], 0 ) ;PRINT_LAST_ERROR;


    PetscPrintf(PETSC_COMM_SELF, "[%D]\tUpdating shared mesh.\n", rank);
    result = parallel_communicators[isystem]->update_shared_mesh() ;PRINT_LAST_ERROR;

    //PetscPrintf(PETSC_COMM_SELF, "[%D]\tResolving shared ents and ghosts.\n", rank);
    //result = parallel_communicators[isystem]->exchange_ghost_cells(dim, 0, 1, 1, true);PRINT_LAST_ERROR;
    //result = parallel_communicators[isystem]->resolve_shared_ents(root_sets[isystem], dim, 0, &globalId);PRINT_LAST_ERROR;

    //PetscPrintf(PETSC_COMM_SELF, "[%D]\tAssigning global ranks.\n", rank);
  	result = parallel_communicators[isystem]->assign_global_ids(root_sets[isystem], dim, 1, false, false, false);PRINT_LAST_ERROR;


    //mdbImpl->print_database() ;

    PetscPrintf(PETSC_COMM_SELF, "[%D]\tChecking adjacencies.\n", rank);
    result = mdbImpl->check_adjacencies() ;PRINT_LAST_ERROR;
  }

	PetscBarrier(0);

  if (build_solutions)
  {
    // Build the nodal solution values & get the variable
    // names from the EquationSystems object
    std::vector<Number>      soln;
    std::vector<std::string> names;

    PetscPrintf(PETSC_COMM_SELF, "[%D]\tBuilding solution.\n", rank);
    const unsigned int n_es_vars = physics.n_vars();
    physics.build_variable_names(names);

    physics.build_solution_vector(soln);

    if (!rank)
    {
      std::vector<Tag> solution_tags(n_es_vars);
      std::vector<double> var_solution_nodal(n_nodes);

      for (unsigned int ivars = 0 ; ivars < n_es_vars; ivars++)
      {
        result = mdbImpl->tag_get_handle( names[ivars].c_str(), 1, MB_TYPE_DOUBLE,
                          solution_tags[ivars], MB_TAG_DENSE|MB_TAG_CREAT );PRINT_LAST_ERROR;

        for(node_it = mesh.nodes_begin(), index=0;node_it!=node_end;++node_it, ++index)
        {
          unsigned int nd_id = (*node_it)->id();

          var_solution_nodal[index] = soln[nd_id*n_es_vars+ivars] ;

          // PetscPrintf (PETSC_COMM_WORLD,"\n Variable %s\tNode %d\tValue %10.7e ", names[ivars].c_str(), nd_id, var_solution_nodal[index]);
        }

        result = mdbImpl->tag_set_data(solution_tags[ivars], &global_nodes[0], n_nodes, &var_solution_nodal[0]);PRINT_LAST_ERROR;
      }

      if (verbose)
      {
        PetscPrintf(PETSC_COMM_WORLD, "[%D]\tWriting file.\n", rank);

        iMesh_save(imesh, (iBase_EntitySetHandle)root_sets[isystem], fname.c_str(), writeOpts.c_str(),
                         (int*)(&result), fname.length(), writeOpts.length()); PRINT_LAST_ERROR;

        PetscPrintf(PETSC_COMM_WORLD, "[%D]\tWriting file done.\n", rank);
      }

    }

  	PetscBarrier(0);
  }


	return result;

}



ErrorCode check_consistent_ids( Interface* mdbImpl,
                                const EntityHandle* entities,
                                const int* orig_ids,
                                int num_ents,
                                const char* singular_name,
                                const char* plural_name )
{
  ErrorCode result;
  int rank, size, ierr;
  MPI_Comm_rank( PETSC_COMM_WORLD, &rank );
  MPI_Comm_size( PETSC_COMM_WORLD, &size );

  Tag id_tag;
  result = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, id_tag ); PRINT_LAST_ERROR;
  std::vector<int> new_ids(num_ents);
  result = mdbImpl->tag_get_data( id_tag, entities, num_ents, &new_ids[0] ); PRINT_LAST_ERROR;
  // This test is wrong.  a) The caller can select a start ID so there's
  // no guarantee that the IDs will be in any specific range and b) There
  // is no reason to expect the global number of entities to be num_ents*size
  // if there are any shared entities. J.Kraftcheck 2010-12-22
  //rval = MB_SUCCESS;
  //for (int i = 0; i < num_ents; ++i)
  //  if (new_ids[i] < 0 || new_ids[i] >= num_ents*size) {
  //    std::cerr << "ID out of bounds on proc " << rank
  //              << " : " << new_ids[i] << " not in [0," << num_ents*size-1
  //              << "]" << std::endl;
  //    rval = MB_FAILURE;
  //  }
  //if (MB_SUCCESS != rval)
  //  return rval;

    // Gather up all data on root proc for consistency check
  std::vector<int> all_orig_ids(num_ents*size), all_new_ids(num_ents*size);
  ierr = MPI_Gather( (void*)orig_ids, num_ents, MPI_INT, &all_orig_ids[0], num_ents, MPI_INT, 0, PETSC_COMM_WORLD );
  if (ierr)
    return MB_FAILURE;
  ierr = MPI_Gather( &new_ids[0], num_ents, MPI_INT, &all_new_ids[0], num_ents, MPI_INT, 0, PETSC_COMM_WORLD );
  if (ierr)
    return MB_FAILURE;

    // build a local map from original ID to new ID and use it
    // to check for consistancy between all procs
  result = MB_SUCCESS;;
  if (0 == rank) {
      // check for two processors having different global ID for same entity
    std::map<int,int> idmap; // index by original ID and contains new ID
    std::map<int,int> owner; // index by original ID and contains owning rank
    for (int i = 0; i < num_ents*size; ++i) {
      std::map<int,int>::iterator it = idmap.find(all_orig_ids[i]);
      if (it == idmap.end()) {
        idmap[all_orig_ids[i]] = all_new_ids[i];
        owner[all_orig_ids[i]] = i/num_ents;
        PetscPrintf(PETSC_COMM_WORLD, "\nConsistant %s IDs between processors %D and %D : %D and %D respectively.", singular_name, owner[all_orig_ids[i]], i/num_ents, it->second, all_new_ids[i]);
      }
      else if (it->second != all_new_ids[i]) {
        PetscPrintf(PETSC_COMM_WORLD, "\nInconsistant %s IDs between processors %D and %D : %D and %D respectively.", singular_name, owner[all_orig_ids[i]], i/num_ents, it->second, all_new_ids[i]);
        result = MB_FAILURE;
      }
    }
      // check for two processors having same global ID for different entities
    idmap.clear();
    owner.clear();
    for (int i = 0; i < num_ents*size; ++i) {
      std::map<int,int>::iterator it = idmap.find(all_new_ids[i]);
      if (it == idmap.end()) {
        idmap[all_new_ids[i]] = all_orig_ids[i];
        owner[all_new_ids[i]] = i/num_ents;
      }
      else if (it->second != all_orig_ids[i]) {
        std::cerr << "ID " << all_new_ids[i]
                  << " assigned to different " << plural_name << " on processors "
                  << owner[all_new_ids[i]] << " and " << i/num_ents << std::endl;
        result = MB_FAILURE;
      }
    }
  }
  return result;
}


ErrorCode MOAB_libMesh_Wrapper::write_equation_systems_parallel(const std::string& fname, const EquationSystems& physics, unsigned int isystem,
                            bool build_solutions, bool build_adjacencies, bool verbose)
{
	ErrorCode result = MB_SUCCESS;
	int nprocs, rank;
	std::string writeOpts = "PARALLEL=WRITE_PART;DEBUG_IO=0" ;

	assert(isystem < n_systems);

	MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

	const MeshBase& mesh = (MeshBase&)physics.get_mesh();
	unsigned int dim = mesh.spatial_dimension();

	PerfLog logger("MOAB_libMesh_Wrapper");

	logger.push("bootstrap()");

	// Initialize the moab object
	result = mdbImpl->set_dimension(dim);

	unsigned int node_offset = 0;
	unsigned int elem_offset = 0;
	for (int iproc=0; iproc < rank; ++iproc)
	{
		node_offset += mesh.n_nodes_on_proc(iproc);
		elem_offset += mesh.n_elem_on_proc(iproc);
	}

	const unsigned int n_nodes = mesh.n_nodes();
	const unsigned int n_local_nodes = mesh.n_local_nodes();
	const unsigned int n_elems = mesh.n_elem();
	const unsigned int n_local_elems = mesh.n_local_elem();
	unsigned int index = 0;

	std::vector<EntityHandle> global_nodes, local_nodes;
	std::vector<EntityHandle> tglobal_nodes;
	std::vector<int> nd_global_ids;

	tglobal_nodes.assign(n_nodes, n_local_nodes);
	global_nodes.assign(n_nodes, 0);
	nd_global_ids.assign(n_local_nodes, 0);

	logger.pop("bootstrap()");

	logger.push("get_nodes()");

  // write mesh nodes onto MOAB instance
  double vcoords[3];
  for(MeshBase::const_node_iterator node_it = mesh.local_nodes_begin(); node_it!=mesh.local_nodes_end(); ++node_it)
  {
    const Node& nd = *(*node_it);

    // store the vertex corresponding to a certain partition
    if (nd.processor_id() == rank)
    {
      vcoords[0] = nd(0); vcoords[1] = nd(1); vcoords[2] = nd(2);
      result = mdbImpl->create_vertex(vcoords, tglobal_nodes[nd.id()]); PRINT_LAST_ERROR;

      local_nodes.push_back(tglobal_nodes[nd.id()]);
      nd_global_ids.push_back(tglobal_nodes[nd.id()]+node_offset);
    }
  }

  assert(local_nodes.size() == n_local_nodes);
  //result = mdbImpl->add_entities ( root_sets[isystem], &local_nodes[0], local_nodes.size() ); PRINT_LAST_ERROR;
  result = mdbImpl->add_entities ( partition_handles[isystem], &local_nodes[0], local_nodes.size() ); PRINT_LAST_ERROR;

  PetscBarrier(0);

//	int ierr = MPI_Allgather(&local_nodes[0], local_nodes.size(), MPI_DOUBLE,
//	              &global_nodes[0], global_nodes.size(), MPI_DOUBLE, PETSC_COMM_WORLD);
  int ierr = MPI_Allreduce(&tglobal_nodes[0], &global_nodes[0], global_nodes.size(), MPI_DOUBLE,
              MPI_SUM, PETSC_COMM_WORLD);
  CHKERRABORT(PETSC_COMM_WORLD,ierr);

  logger.pop("get_nodes()");

  if (verbose && true)
  {
    if (rank)
    {
      for(unsigned int i=0; i < global_nodes.size();++i)
        PetscPrintf(PETSC_COMM_SELF, "\n%D Owned=%D \t tnodes_vtk_id[%D] = %D\t nodes_vtk_id[%D] = %D", rank, mesh.node(i).processor_id()==rank, i, tglobal_nodes[i], i, global_nodes[i]);
    }

    PetscBarrier(0);

    if (!rank)
    {
      for(unsigned int i=0; i < global_nodes.size();++i)
        PetscPrintf(PETSC_COMM_SELF, "\n%D Owned=%D \t tnodes_vtk_id[%D] = %D\t nodes_vtk_id[%D] = %D", rank, mesh.node(i).processor_id()==rank, i, tglobal_nodes[i], i, global_nodes[i]);
    }
  }

	PetscBarrier(0);

	  Range& part_set = parallel_communicators[isystem]->partition_sets();

  logger.push("get_elements()");


  std::vector<unsigned int> cell_connectivity;
  std::vector<EntityHandle> moab_connectivity;

  std::vector<EntityHandle> local_elems;
  std::vector<int> el_global_ids;

  for(MeshBase::const_element_iterator element_it = mesh.active_local_elements_begin(); element_it!=mesh.active_local_elements_end(); ++element_it)
  {
    const Elem* elem = *element_it;

      this->cell_connectivity(elem, cell_connectivity);
      moab_connectivity.resize(cell_connectivity.size());

      for (unsigned int i=0; i < cell_connectivity.size(); ++i) {
        moab_connectivity[i] = tglobal_nodes[cell_connectivity[i]];
      }

      EntityHandle elem_handle ;
      result = mdbImpl->create_element (cell_type(elem), &moab_connectivity[0], moab_connectivity.size(), elem_handle);PRINT_LAST_ERROR;

      local_elems.push_back(elem_handle);
      el_global_ids.push_back(elem_handle);

  //    part_set.insert(elem_handle);
  }
  //result = mdbImpl->add_entities ( root_sets[isystem], &local_elems[0], local_elems.size() ); PRINT_LAST_ERROR;
  result = mdbImpl->add_entities ( partition_handles[isystem], &local_elems[0], local_elems.size() ); PRINT_LAST_ERROR;

  Range ents;
  result = mdbImpl->get_entities_by_type(0, MBENTITYSET, ents, false);
  PetscPrintf(PETSC_COMM_SELF, "\n[%D]\t Number of EntitySets: %D\n", rank, ents.size());

  logger.pop("get_elements()");

  logger.push("set_global_id()");

  // next get the global ID tag handle
  Tag globalId = mdbImpl->globalId_tag();
  result = mdbImpl->tag_set_data(globalId, &local_nodes[0], local_nodes.size(), &nd_global_ids[0]);PRINT_LAST_ERROR;
  result = mdbImpl->tag_set_data(globalId, &local_elems[0], local_elems.size(), &el_global_ids[0]);PRINT_LAST_ERROR;

//    result = check_consistent_ids(mdbImpl, &local_nodes[0], &nd_global_ids[0], local_nodes.size(), "vertex", "vertices");PRINT_LAST_ERROR;
//    result = check_consistent_ids(mdbImpl, &local_elems[0], &el_global_ids[0], local_elems.size(), "element", "elements");PRINT_LAST_ERROR;

  //result = parallel_communicators[isystem]->exchange_ghost_cells(dim, 0, 1, 1, true);PRINT_LAST_ERROR;

  logger.pop("set_global_id()");

  logger.push("make_sets()");

  // create some sets
  EntityHandle vertex_set, region_set, face_set;
  vertex_set = make_set( MESHSET_ORDERED|MESHSET_TRACK_OWNER, &local_nodes[0], local_nodes.size(), false,  MOAB_libMesh_Wrapper::VERTEX_SET_ID*(isystem+1)  );
  region_set = make_set( MESHSET_SET|MESHSET_TRACK_OWNER,     &local_elems[0], local_elems.size(), false, MOAB_libMesh_Wrapper::REGION_SET_ID*(isystem+1)  );

  logger.pop("make_sets()");

  logger.push("set_geom_dimension()");

  // next get the global ID tag handle
  Tag dimId = mdbImpl->geom_dimension_tag();

  int vdim = 0;
  result = mdbImpl->tag_set_data(dimId, &region_set, 1, &dim);PRINT_LAST_ERROR;
  result = mdbImpl->tag_set_data(dimId, &vertex_set, 1, &vdim);PRINT_LAST_ERROR;

  logger.pop("set_geom_dimension()");

  part_set.insert(region_set);
  //part_set.insert(vertex_set);

  if (build_adjacencies)
  {
	logger.push("build_adjacencies()");

    Range dim3, dim2, adj;
    PetscPrintf(PETSC_COMM_SELF, "[%D]\tGetting adjacencies.\n", rank);
    //result = parallel_communicators[isystem]->get_part_entities(dim3, 3);PRINT_LAST_ERROR;
    result = mdbImpl->get_entities_by_dimension( region_set, 3, dim3, true );PRINT_LAST_ERROR;
    PetscPrintf(PETSC_COMM_SELF, "[%D]\t\t%d %d-Dim quantities.\n", rank, dim3.size(), 3);

    // generate all lower dimensional adjacencies
    //result = parallel_communicators[isystem]->get_part_entities(dim2, 2);PRINT_LAST_ERROR;
    result = mdbImpl->get_entities_by_dimension( region_set, 2, dim2, true );PRINT_LAST_ERROR;
    PetscPrintf(PETSC_COMM_SELF, "[%D]\t\t%d %d-Dim quantities.\n", rank, dim2.size(), 2);
    mdbImpl->get_adjacencies( dim3, 1, true, adj, Interface::UNION );
    mdbImpl->get_adjacencies( dim2, 1, true, adj, Interface::UNION );
    mdbImpl->get_adjacencies( dim3, 2, true, adj, Interface::UNION );
    PetscPrintf(PETSC_COMM_SELF, "[%D]\t\t%d adjacency quantities.\n", rank, adj.size());

    mdbImpl->add_entities( partition_handles[isystem], adj );

    int fdim=dim-1;
    face_set   = make_set( MESHSET_SET|MESHSET_TRACK_OWNER,     adj/*dim2*/, false, MOAB_libMesh_Wrapper::FACE_SET_ID*(isystem+1)    );
    //result = mdbImpl->add_entities( partition_handles[isystem], &face_set, 1 ) ;PRINT_LAST_ERROR;
    result = mdbImpl->tag_set_data(dimId, &face_set, 1, &fdim);PRINT_LAST_ERROR;
    //part_set.insert(face_set);

	logger.pop("build_adjacencies()");
  }

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tTotal number of Partition sets %D.\n", rank, part_set.size());

	logger.push("set_miscellaneous_tags()");
//  result = mdbImpl->add_parent_meshset( vertex_set, root_sets[isystem] ) ;PRINT_LAST_ERROR;
//  result = mdbImpl->add_parent_meshset( region_set, root_sets[isystem] ) ;PRINT_LAST_ERROR;

  result = mdbImpl->add_parent_meshset( partition_handles[isystem], root_sets[isystem] ) ;PRINT_LAST_ERROR;
  //result = mdbImpl->add_entities( root_sets[isystem], &partition_handles[isystem], 1 ) ;PRINT_LAST_ERROR;
  result = mdbImpl->unite_meshset ( root_sets[isystem], partition_handles[isystem] ); PRINT_LAST_ERROR;

  Tag partn_tag = parallel_communicators[isystem]->partition_tag	();
  // set the partition to the global vertex/region set => 1 proc file
  result = mdbImpl->tag_set_data(partn_tag, &partition_handles[isystem], 1, &rank);PRINT_LAST_ERROR;

  result = parallel_communicators[isystem]->set_partitioning(partition_handles[isystem]);PRINT_LAST_ERROR;

  logger.pop("set_miscellaneous_tags()");

  logger.push("assign_global_ids()");
  PetscPrintf(PETSC_COMM_SELF, "[%D]\tAssigning global ranks.\n", rank);
  result = parallel_communicators[isystem]->assign_global_ids(root_sets[isystem], dim, 1, false, false, false);PRINT_LAST_ERROR;
  logger.pop("assign_global_ids()");

  logger.push("resolve_shared_ents()");
  PetscPrintf(PETSC_COMM_SELF, "[%D]\tResolving shared ents and ghosts. Partition Handle %D\n", rank, partition_handles[isystem]);
  result = parallel_communicators[isystem]->resolve_shared_sets(root_sets[isystem]);PRINT_LAST_ERROR;
  result = parallel_communicators[isystem]->resolve_shared_ents(root_sets[isystem], dim, 0 /*, &globalId*/);PRINT_LAST_ERROR;
  result = parallel_communicators[isystem]->exchange_ghost_cells(dim, 0, 1, 1, false, true, &root_sets[isystem]);PRINT_LAST_ERROR;
  logger.pop("resolve_shared_ents()");

  logger.push("check_all_shared_handles()");
  result = parallel_communicators[isystem]->check_all_shared_handles();PRINT_LAST_ERROR;
  logger.pop("check_all_shared_handles()");

  if (verbose) {
    moab::Range ents;
    // iterate over dimensions
    for (unsigned int d = 0; d <= dim; d++) {
      ents.clear();
      result = mdbImpl->get_entities_by_dimension((EntityHandle)root_sets[isystem], d, ents);PRINT_LAST_ERROR;
      PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tFound %d %d-dimensional entities\n", rank, ents.size(), d);
      /*
      for (moab::Range::iterator it = ents.begin(); it != ents.end(); it++) {
        moab::EntityHandle ent = *it;
        PetscPrintf(PETSC_COMM_SELF, "[%D]\t\tFound d=%d entity %d.\n", rank, d, mdbImpl->id_from_handle(ent));
      }
      */
    }
  }

//    result = mdbImpl->add_entities( root_sets[isystem], &region_set, 1 ) ;PRINT_LAST_ERROR;
//  	result = mdbImpl->add_entities( root_sets[isystem], &vertex_set, 1 ) ;PRINT_LAST_ERROR;

  // finally add the file set to the mesh's root set
  //result = mdbImpl->add_entities( 0, &root_sets[isystem], 1 ) ;PRINT_LAST_ERROR;
  //result = mdbImpl->add_parent_meshset( root_sets[isystem], 0 ) ;PRINT_LAST_ERROR;


  PetscPrintf(PETSC_COMM_SELF, "[%D]\tUpdating shared mesh.\n", rank);
  logger.push("update_shared_mesh()");
  //result = parallel_communicators[isystem]->create_interface_sets(dim, dim-1) ;PRINT_LAST_ERROR;
  result = parallel_communicators[isystem]->update_shared_mesh() ;PRINT_LAST_ERROR;
  result = parallel_communicators[isystem]->collective_sync_partition() ;PRINT_LAST_ERROR;
  logger.pop("update_shared_mesh()");

  if (verbose && false)
  {
    if (!rank)
      mdbImpl->print_database() ;

    PetscBarrier(0);

    if (rank)
      mdbImpl->print_database() ;

	  PetscBarrier(0);
	}

  if (build_adjacencies && verbose)
  {
	logger.push("check_adjacencies()");
    PetscPrintf(PETSC_COMM_SELF, "[%D]\tChecking adjacencies.\n", rank);
    result = mdbImpl->check_adjacencies() ;PRINT_LAST_ERROR;
	logger.pop("check_adjacencies()");
  }


	PetscBarrier(0);

  if (build_solutions)
  {
	logger.push("build_solutions()");
    // Build the nodal solution values & get the variable
    // names from the EquationSystems object
    std::vector<Number>      soln;
    std::vector<std::string> names;

    PetscPrintf(PETSC_COMM_SELF, "[%D]\tBuilding solution.\n", rank);
    const unsigned int n_es_vars = physics.n_vars();
    physics.build_variable_names(names);

    physics.build_solution_vector(soln);

    std::vector<Tag> solution_tags(n_es_vars);
    std::vector<double> var_solution_nodal(n_nodes);

    for (unsigned int ivars = 0 ; ivars < n_es_vars; ivars++)
    {
      result = mdbImpl->tag_get_handle( names[ivars].c_str(), 1, MB_TYPE_DOUBLE,
                        solution_tags[ivars], MB_TAG_DENSE|MB_TAG_CREAT );PRINT_LAST_ERROR;

      for(MeshBase::const_node_iterator node_it = mesh.local_nodes_begin(); node_it!=mesh.local_nodes_end(); ++node_it)
      //for(node_it = mesh.local_nodes_begin(), index=0;node_it!=node_end;++node_it, ++index)
      {
        unsigned int nd_id = (*node_it)->id();

        var_solution_nodal[index] = soln[nd_id*n_es_vars+ivars] ;

        // PetscPrintf (PETSC_COMM_WORLD,"\n Variable %s\tNode %d\tValue %10.7e ", names[ivars].c_str(), nd_id, var_solution_nodal[index]);
      }

      result = mdbImpl->tag_set_data(solution_tags[ivars], &local_nodes[0], n_local_nodes, &var_solution_nodal[0]);PRINT_LAST_ERROR;
    }
	logger.pop("build_solutions()");

    if (verbose)
    {
		logger.push("write_file()");
      PetscPrintf(PETSC_COMM_WORLD, "[%D]\tWriting file.\n", rank);

      result = mdbImpl->write_file(fname.c_str(), "MOAB", writeOpts.c_str());PRINT_LAST_ERROR;
//      iMesh_save(imesh, (iBase_EntitySetHandle)root_sets[isystem], fname.c_str(), writeOpts.c_str(),
//                       (int*)(&result), fname.length(), writeOpts.length()); PRINT_LAST_ERROR;

      PetscPrintf(PETSC_COMM_WORLD, "[%D]\tWriting file done.\n", rank);
		logger.pop("write_file()");
    }

  	PetscBarrier(0);
  }

	return result ;

}



ErrorCode MOAB_libMesh_Wrapper::write_equation_systems_parallel_old(const std::string& fname, const EquationSystems& physics, unsigned int isystem)
{
	ErrorCode result;
	int nprocs, rank;

	MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

	//std::string par_writeOpts(" moab:PARALLEL=WRITE_PART moab:CPUTIME ");
	std::string par_writeOpts(" moab:PARALLEL=WRITE_PART moab:CPUTIME ");
	std::string ser_writeOpts("");
	std::string writeOpts = (nprocs > 1 ? par_writeOpts : ser_writeOpts) ;

  parallel_communicators[isystem] = moab::ParallelComm::get_pcomm(mdbImpl, partition_handles[isystem], &PETSC_COMM_WORLD);

	ParallelComm& pcomm = *parallel_communicators[isystem];

	const MeshBase& mesh = (MeshBase&)physics.get_mesh();
	unsigned int dim = mesh.spatial_dimension();
	MeshBase::const_element_iterator element_it = mesh.active_local_elements_begin() ;
	const MeshBase::const_element_iterator element_end = mesh.active_local_elements_end();
	MeshBase::const_node_iterator node_it = mesh.local_nodes_begin() ;
	const MeshBase::const_node_iterator node_end = mesh.nodes_end();
	const MeshBase::const_node_iterator node_lend = mesh.local_nodes_end();

	parallel_only() ;

	const unsigned int n_nodes = mesh.parallel_n_nodes();
	const unsigned int n_local_nodes = mesh.n_local_nodes();
	const unsigned int n_elems = mesh.parallel_n_elem();
	const unsigned int n_local_elems = mesh.n_local_elem();
	unsigned int index = 0;

	unsigned int node_offset = 0;
	for (int iproc=0; iproc < rank; ++iproc)
		node_offset += mesh.n_nodes_on_proc(iproc);

	PetscPrintf(PETSC_COMM_SELF, "\n[%D] Node offset = %D", rank, node_offset);

	std::vector<EntityHandle> nodes_vtk_id, tnodes_vtk_id, local_nodes_vtk_id;
	tnodes_vtk_id.assign(n_nodes, 0);
//	local_nodes_vtk_id.assign(n_local_nodes, 0);

	result = mdbImpl->create_meshset( MESHSET_SET|MESHSET_TRACK_OWNER, root_sets[isystem] ); PRINT_LAST_ERROR;
	//root_sets[isystem] = make_set( MESHSET_SET, 0, 0, true,  EMPTY_SET_ID  );

	// Initialize the moab object
	result = mdbImpl->set_dimension(dim);

	PetscBarrier(0);

	// next get the global ID tag handle
	Tag dimId = mdbImpl->geom_dimension_tag();

	// write mesh nodes onto MOAB instance
	double vcoords[3];
//	for(node_it = mesh.local_nodes_begin(), index=0; node_it!=node_lend; ++node_it,++index)
	for(node_it = mesh.nodes_begin(), index=0; node_it!=node_end; ++node_it)
	{
	  const Node& nd = *(*node_it);
	  if (nd.processor_id() == rank)
	  {
      vcoords[0] = nd(0); vcoords[1] = nd(1); vcoords[2] = nd(2);
      //PetscPrintf(PETSC_COMM_SELF, "\n[%D] %D \tNode offset = %D, Node ID = %D", rank, index, node_offset, node_offset+nd.id());
      //nodes_vtk_id[index] = (EntityHandle)(nd.id()) ;
      result = mdbImpl->create_vertex(vcoords, tnodes_vtk_id[node_offset+index]); PRINT_LAST_ERROR;

//	  		PetscPrintf(PETSC_COMM_SELF, "\n[%D]\t Node %D", rank, tnodes_vtk_id[node_offset+index]) ;
	  if (nd.processor_id() == rank)
	    local_nodes_vtk_id.push_back(node_offset+tnodes_vtk_id[node_offset+index]);

	    ++index;

//	  result = mdbImpl->add_parent_meshset (nodes_vtk_id[index], root_sets[isystem]); PRINT_LAST_ERROR;
	  }
//	  else
//	    tnodes_vtk_id[index] = index+1;

	}

	PetscBarrier(0);

//	nodes_vtk_id = tnodes_vtk_id;
  nodes_vtk_id.resize(n_nodes);

//	int err = MPI_Allgather(&local_nodes_vtk_id[0], local_nodes_vtk_id.size(), MPI_DOUBLE,
//	              &nodes_vtk_id[0], local_nodes_vtk_id.size(), MPI_DOUBLE, PETSC_COMM_WORLD);
	int err = MPI_Allreduce(&tnodes_vtk_id[0], &nodes_vtk_id[0], nodes_vtk_id.size(), MPI_DOUBLE,
	              MPI_SUM, PETSC_COMM_WORLD); CHKERRABORT(PETSC_COMM_WORLD,err);


if (rank)
{
  for(unsigned int i=0; i < nodes_vtk_id.size();++i)
    PetscPrintf(PETSC_COMM_SELF, "\n%D tnodes_vtk_id[%D] = %D\t nodes_vtk_id[%D] = %D", rank, i, tnodes_vtk_id[i], i, nodes_vtk_id[i]);
  }

	PetscBarrier(0);

if (!rank)
{
  for(unsigned int i=0; i < nodes_vtk_id.size();++i)
    PetscPrintf(PETSC_COMM_SELF, "\n%D tnodes_vtk_id[%D] = %D\t nodes_vtk_id[%D] = %D", rank, i, tnodes_vtk_id[i], i, nodes_vtk_id[i]);
  }

	PetscBarrier(0);


	// next get the global ID tag handle
	Tag globalId = mdbImpl->globalId_tag();
//	result = mdbImpl->tag_set_data(globalId, &nodes_vtk_id[0], n_nodes, &nodes_vtk_id[0]);PRINT_LAST_ERROR;
	result = mdbImpl->tag_set_data(globalId, &nodes_vtk_id[node_offset], local_nodes_vtk_id.size(), &local_nodes_vtk_id[0]);PRINT_LAST_ERROR;
	//result = mdbImpl->tag_set_data(globalId, &local_elems_vtk_id[0], local_elems_vtk_id.size(), &local_elems_vtk_id[0]);PRINT_LAST_ERROR;


//	iMesh_save(imesh, (iBase_EntitySetHandle)root_sets[isystem], "verts.h5m", writeOpts.c_str(),
//                   (int*)(&result), strlen("verts.h5m"), writeOpts.length()); PRINT_LAST_ERROR;


	assert(local_nodes_vtk_id.size() == n_local_nodes);
	result = mdbImpl->add_entities ( root_sets[isystem], &local_nodes_vtk_id[0], local_nodes_vtk_id.size() ); PRINT_LAST_ERROR;

	//std::vector<int> entity_dims;
	//entity_dims.assign ( n_nodes, 0 );
	//result = mdbImpl->tag_set_data(dimId, &nodes_vtk_id[node_offset], n_nodes, &entity_dims);PRINT_LAST_ERROR;


  std::vector<unsigned int> vtk_cell_connectivity;
  std::vector<EntityHandle> moab_connectivity;

	std::vector<EntityHandle> local_elems_vtk_id;

	for(element_it = mesh.elements_begin(), index=0; element_it!=element_end; ++element_it,++index)
//	for(element_it = mesh.active_local_elements_begin(), index=0; element_it!=mesh.active_local_elements_end(); ++element_it,++index)
	{
		const Elem* elem = *element_it;

		if (elem->processor_id() == rank)
		{
      this->cell_connectivity(elem, vtk_cell_connectivity);
      moab_connectivity.resize(vtk_cell_connectivity.size());

      PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tConnectivity Element %D", rank, index) ;
      for (unsigned int i=0; i < vtk_cell_connectivity.size(); ++i) {
        moab_connectivity[i] = nodes_vtk_id[vtk_cell_connectivity[i]];
  //			moab_connectivity[i] = vtk_cell_connectivity[i]+1;
        PetscPrintf(PETSC_COMM_SELF, "\t%D %D", vtk_cell_connectivity[i], moab_connectivity[i]) ;
      }

      EntityHandle elem_handle ;
      result = mdbImpl->create_element (cell_type(elem), &moab_connectivity[0], moab_connectivity.size(), elem_handle);PRINT_LAST_ERROR;

      // add the adjacency information for each node in the element to other nodes
      result = mdbImpl->add_adjacencies(elem_handle, &moab_connectivity[0], moab_connectivity.size(), true);

	    local_elems_vtk_id.push_back(elem_handle);
	  }

	}

	assert(local_elems_vtk_id.size() == n_local_elems);
	result = mdbImpl->add_entities ( root_sets[isystem], &local_elems_vtk_id[0], local_elems_vtk_id.size() ); PRINT_LAST_ERROR;

/*
	for(element_it = mesh.active_local_elements_begin(), index=0; element_it!=mesh.active_local_elements_end(); ++element_it,++index)
	{
		const Elem* elem = *element_it;
		unsigned int nside_index = 0;
		std::vector<EntityHandle> neighbor_ids(elem->n_sides());
		//neighbor_ids.push_back(elems_vtk_id[index]);
		for (unsigned int side=0; side<elem->n_sides(); ++side)
		{
			if (elem->neighbor(side) != NULL)	// we have neighbours
			{
				neighbor_ids[nside_index] = (EntityHandle)(elem->neighbor(side)->id() );

				std::vector<EntityHandle> adjacencies;
				EntityHandle from_entities[2] = {elems_vtk_id[index], neighbor_ids[nside_index]};
				// generate all edges for these two hexes
				//result = mdbImpl->get_adjacencies( from_entities, 2, 0, true, adjacencies, Interface::INTERSECT); PRINT_LAST_ERROR; //Interface::UNION
				//adjacencies.clear();
				// now find the edges common to both hexes
				//result = mdbImpl->get_adjacencies( from_entities, 2, 2, true, adjacencies, Interface::UNION); PRINT_LAST_ERROR;

				nside_index++;

			}
		}
		//result = mdbImpl->add_adjacencies(elems_vtk_id[index], &neighbor_ids[0], nside_index, true);PRINT_LAST_ERROR;
		result = mdbImpl->get_adjacencies(elems_vtk_id[index], &neighbor_ids[0], nside_index, true);PRINT_LAST_ERROR;

	}
	*/

	// Set the dimension tag for all elements and add to the root set
	//result = mdbImpl->add_entities ( root_sets[isystem], &elems_vtk_id[0], elems_vtk_id.size() ); PRINT_LAST_ERROR;
	//std::vector<int> entity_dims;
	//entity_dims.assign ( local_elems_vtk_id.size(), dim );
	//result = mdbImpl->tag_set_data(dimId, &local_elems_vtk_id[0], local_elems_vtk_id.size(), &entity_dims[0]);PRINT_LAST_ERROR;

	// create some sets
	EntityHandle vertex_set, region_set/*, face_set, empty_set*/;

	//face_set   = make_set( MESHSET_SET,     faces,    12, false, FACE_SET_ID    );
	//vertex_set = make_set( MESHSET_ORDERED|MESHSET_TRACK_OWNER, &nodes_vtk_id[node_offset], n_local_nodes, true,  VERTEX_SET_ID  );

			PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tLocal nodes: %D\tLocal elements: %D", rank, local_nodes_vtk_id.size(), local_elems_vtk_id.size()) ;

	vertex_set = make_set( MESHSET_SET|MESHSET_TRACK_OWNER, &local_nodes_vtk_id[0], local_nodes_vtk_id.size(), false,  VERTEX_SET_ID  );
	region_set = make_set( MESHSET_SET|MESHSET_TRACK_OWNER,     &local_elems_vtk_id[0], local_elems_vtk_id.size(), false, REGION_SET_ID  );

//	vertex_set = make_set( MESHSET_ORDERED|MESHSET_TRACK_OWNER, &nodes_vtk_id[0], nodes_vtk_id.size(), false,  VERTEX_SET_ID  );
//	region_set = make_set( MESHSET_SET|MESHSET_TRACK_OWNER,     &elems_vtk_id[0], elems_vtk_id.size(), false, REGION_SET_ID  );

	//const unsigned empty_flags = MESHSET_ORDERED|MESHSET_TRACK_OWNER;
	//empty_set  = make_set( empty_flags,     0,         0, true,  EMPTY_SET_ID   );
	//EntityHandle sets[] = {face_set, vertex_set, region_set, empty_set};
	//EntityHandle sets[] = { vertex_set, region_set, empty_set};
	//make_set( MESHSET_ORDERED, sets,      3, false, SET_SET_ID     );

	PetscBarrier(0);

//	if (rank) {
	// ! element dimensions seems to be screwed up !
	//result = mdbImpl->create_meshset( MESHSET_SET, region_set );PRINT_LAST_ERROR;
	//for (unsigned int i = 0; i < n_local_nodes; ++i)
	//	PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tNode id = %D", rank, nodes_vtk_id[node_offset+i]) ;
	//for (unsigned int i = 0; i < n_local_elems; ++i) {
		//PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tElement id = %D", rank, elems_vtk_id[i]) ;
		//mdbImpl->list_entity(elems_vtk_id[i]);
		//// result = mdbImpl->add_parent_child( region_set, elems_vtk_id[i] );PRINT_LAST_ERROR;
	//}
	//}
	//PetscBarrier(0);


  bool verbose = false;
  if (verbose) {
  	PetscBarrier(0);

    moab::Range ents;
    // iterate over dimensions
    for (unsigned int d = 0; d <= dim; d++) {
      ents.clear();
      result = mdbImpl->get_entities_by_dimension((EntityHandle)root_sets[isystem], d, ents);PRINT_LAST_ERROR;
      PetscPrintf(PETSC_COMM_SELF, "\n[%D]\tFound %d %d-dimensional entities\n", rank, ents.size(), d);
      for (moab::Range::iterator it = ents.begin(); it != ents.end(); it++) {
        moab::EntityHandle ent = *it;
        PetscPrintf(PETSC_COMM_SELF, "[%D]\t\tFound d=%d entity %d.\n", rank, d, mdbImpl->id_from_handle(ent));
      }
    }
	}


	// create some set parent-child links
	//result = iface->add_parent_child( face_set, vertex_set) ;PRINT_LAST_ERROR;
	//result = iface->add_child_meshset( region_set, face_set ) ;PRINT_LAST_ERROR;

//	result = mdbImpl->add_parent_meshset( vertex_set, region_set ) ;PRINT_LAST_ERROR;
//	result = mdbImpl->add_parent_meshset( region_set, root_sets[isystem] ) ;PRINT_LAST_ERROR;
//	result = mdbImpl->add_parent_meshset( vertex_set, root_sets[isystem] ) ;PRINT_LAST_ERROR;

	PetscBarrier(0);

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tGetting adjacencies.\n", rank);
	std::vector<EntityHandle> adjacencies2, adjacencies1, adjacencies0;
	//result = mdbImpl->get_adjacencies(&local_elems_vtk_id[0], local_elems_vtk_id.size(), 0, true, adjacencies2, Interface::INTERSECT);PRINT_LAST_ERROR;
	//result = mdbImpl->get_adjacencies(&adjacencies2[0], adjacencies2.size(), 0, true, adjacencies1, Interface::INTERSECT);PRINT_LAST_ERROR;
	//result = mdbImpl->get_adjacencies(&adjacencies1[0], adjacencies1.size(), 0, true, adjacencies0, Interface::INTERSECT);PRINT_LAST_ERROR;


  Tag partn_tag = pcomm.partition_tag	();
	result = mdbImpl->tag_set_data(partn_tag, &region_set, 1, &rank);PRINT_LAST_ERROR;
	result = mdbImpl->tag_set_data(partn_tag, &vertex_set, 1, &rank);PRINT_LAST_ERROR;
	//std::vector<int> partitions(n_local_nodes, rank);
	//result = mdbImpl->tag_set_data(partn_tag, &local_elems_vtk_id[0], local_elems_vtk_id.size(), &partitions[0]);PRINT_LAST_ERROR;
	//result = mdbImpl->tag_set_data(partn_tag, &local_nodes_vtk_id[0], local_nodes_vtk_id.size(), &partitions[0]);PRINT_LAST_ERROR;

	EntityHandle local_vpart, local_epart;
	result = pcomm.set_partitioning(region_set);PRINT_LAST_ERROR;
	result = pcomm.create_part(local_epart);
	result = pcomm.set_partitioning(vertex_set);PRINT_LAST_ERROR;
	result = pcomm.create_part(local_vpart);

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tUpdating shared mesh.\n", rank);
	result = pcomm.update_shared_mesh() ;PRINT_LAST_ERROR;

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tResolving shared ents and ghosts.\n", rank);
	result = pcomm.exchange_ghost_cells(dim, 0, 1, 1, true);PRINT_LAST_ERROR;
	result = pcomm.resolve_shared_ents(root_sets[isystem], dim, 0, &globalId);PRINT_LAST_ERROR;

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tAssigning global ranks.\n", rank);
//	result = pcomm.assign_global_ids(root_sets[isystem], dim, 1, false, true, false);PRINT_LAST_ERROR;

	result = mdbImpl->tag_set_data(dimId, &region_set, 1, &dim);PRINT_LAST_ERROR;
	int vdim = 0;
	result = mdbImpl->tag_set_data(dimId, &vertex_set, 1, &vdim);PRINT_LAST_ERROR;
	//entity_dims.resize( n_local_elems ) ;
	//entity_dims.assign ( n_local_elems, dim );
	//result = mdbImpl->tag_set_data(dimId, &elems_vtk_id[0], n_local_elems, &entity_dims);PRINT_LAST_ERROR;

	if (rank)
		mdbImpl->print_database() ;

	PetscBarrier(0);

if (!rank)
		mdbImpl->print_database() ;

	PetscBarrier(0);

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tChecking adjacencies.\n", rank);
	result = mdbImpl->check_adjacencies() ;PRINT_LAST_ERROR;

  bool build_tags = false;

  if (build_tags)
  {
    // Build the nodal solution values & get the variable
    // names from the EquationSystems object
    std::vector<Number>      soln;
    std::vector<std::string> names;

    PetscPrintf(PETSC_COMM_SELF, "[%D]\tBuilding solution.\n", rank);
    const unsigned int n_es_vars = physics.n_vars();
    physics.build_variable_names(names);

    std::vector<Tag> solution_tags(n_es_vars);
    std::vector<double> var_solution_nodal(n_local_nodes);

    physics.build_solution_vector(soln);

    for (unsigned int ivars = 0 ; ivars < n_es_vars; ivars++)
    {
      result = mdbImpl->tag_get_handle( names[ivars].c_str(), 1, MB_TYPE_DOUBLE,
                        solution_tags[ivars], MB_TAG_DENSE|MB_TAG_CREAT );PRINT_LAST_ERROR;

      for(node_it = mesh.local_nodes_begin(), index=0;node_it!=node_lend;++node_it, ++index)
      {
        unsigned int nd_id = (*node_it)->id();

        var_solution_nodal[index] = soln[nd_id*n_es_vars+ivars] ;

        // fprintf (pFile,"%10.7e ",(*sys.solution)(dof_nr));
      }

      result = mdbImpl->tag_set_data(solution_tags[ivars], &local_nodes_vtk_id[0], n_local_nodes, &var_solution_nodal[0]);PRINT_LAST_ERROR;
    }

    PetscBarrier(0);
	}

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tWriting file.\n", rank);
//	if (!rank)
	{
	iMesh_save(imesh, (iBase_EntitySetHandle)root_sets[isystem], fname.c_str(), writeOpts.c_str(),
                   (int*)(&result), fname.length(), writeOpts.length()); PRINT_LAST_ERROR;
	  //result = mdbImpl->write_file(fname.c_str(), 0, writeOpts.c_str());PRINT_LAST_ERROR;
	}

  PetscPrintf(PETSC_COMM_SELF, "[%D]\tWriting file done.\n", rank);
//	PetscBarrier(0);

	return result;

}
